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The  continuity  equations  of  the  drift-diffusion  semiconductor  device  model  are  hard  to  dis¬ 
cretize  because  of  the  severe  variations  of  variables.  The  Scharfetter  -  Gummel  scheme  applies 
harmonic  averages  to  the  conductance  function,  and  usually  produces  acceptable  solutions  despite 
the  sharp  parameter  changes.  Some  authors  attribute  the  success  of  this  scheme  to  the  relative 
smoothness  of  the  currents,  as  compared  with  carrier  concentrations.  In  this  report  it  is  shown  that 
current  smoothness  cannot  be  derived  from  the  differential  equations,  but  is  related  to  the  specific 
boundary  condition  configuration.  The  success  of  the  Scherfetter  -  Gummel  method  is  hence  easy  to 
justify  for  nearly  1-D  devices,  but  is  harder  to  justify  for  some  other  geometries.  A  stream  function 
formulation  related  to  that  scheme  is  shown  to  overcome  cases  where  the  continuity  equations  are 
ill  -  conditioned. 
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1  introduction 

The  drift-diffusion  semiconductor  equations  are  usually  discretized  by  the 
Scharfetter  -  Gummel  scheme.  This  method  applies  harmonic  averages  to 
the  carrier  conductance  function  [2]. 

Assuming  1-dimensional,  constant  current,  and  piecewise  linear  potential 
variation,  this  scheme  is  exact.  Mock  [6]  shows  that  the  Scharfetter  -  Gum¬ 
mel  discretization  is  equivalent  to  a  certain  finite  -  element  formulation  for  a 
stream  potential.  He  attributes  the  success  of  this  scheme  to  the  experimen¬ 
tal  fact  that  currents  in  semiconductors  are  usually  smoother  than  carrier 
concentrations. 

Brezzi  [3]  shows  the  relaion  of  this  scheme  to  a  hybrid  finite  element  one, 
where  the  currents  are  discretized  by  constant  elements,  and  carrier  con¬ 
centrations  by  linear  ones.  Brezzi  further  states  that  “...  a  solid  theoreti¬ 
cal  argument  showing  cne  superiority  of  the  harmonic  average  is  still  to  be 
found,  in  our  opinion,  except  for  the  obvious  1-D  case..”  In  this  report  it 
is  shown  that  in  the  2-D  case,  the  equation  for  the  electron  stream  function 
is  equivalent  to  the  equation  for  the  hole  Slotbloom  variable,  and  similarly: 
the  equation  for  the  hole  stream  function  is  equivalent  to  the  equation  for 
electron  Slotbloom  variable. 

Smoothness  of  the  stream  function  can  hence  be  attributed  only  to  the  dif¬ 
ferent  boundary  conditions.  The  advantage  of  harmonic  averages  is  conse- 
quantly  easy  to  justify  for  long,  narrow  devices.  It  may  however  not  exsist 
in  some  other  geometrical  configurations. 

While  the  Scherfetter  -  Gummel  discretization  and  the  stream  function 
scheme  produce  the  same  solution  (assuming  infinite  computational  preci¬ 
sion),  it  is  shown  in  section  3  that  the  latter  may  be  far  better  numerically 
conditioned. 

2  Model  and  Discretization 

The  scaled  steady  state  drift  -  diffusion  semiconductor  equations,  expressed 
in  Slotbloom  Variables  are: 

A2At/>  -  e^u  -  e~wv  -  C(x)  =  0  (1) 

V-e0Vu  =  f?  (2) 
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Ve~^Vv  =  R  (3) 

Here  ip  is  scaled  electrostatic  potential,  u  and  v  are  exponentials  of  the 
pseudo-Fermi  potentials,  C  is  the  doping  concentration,  and  R  is  the  recom¬ 
bination  rate.  The  electron  and  hole  concentrations  n,p  are  related  to  u,v 
by  the  relations:  n  =  e^u,  v  =  e~^v. 

For  this  discussion  it  will  be  assumed  that  the  boundary  conditions  for  p,  u,  v 
are  of  Dirichlet  type  on  a  portion  Tp  of  the  boundary  T,  and  of  Neuman 
type  on  T yv  where  T/v  =  F\rx? -  The  Dirichlet  portion  corresponds  to  Ohmic 
boundary  conditions,  and  the  Neumann  portion  to  reflecting  boundaries.  In 
some  cases  the  carrier  concentration  variables  have  internal  Neumann  bound¬ 
aries  which  are  not  on  part  of  Tpj.  This  occurrs  at  insulator  -  semiconductor 
interfaces. 

The  discretization  of  (  2)(  3)  may  introduce  a  large  truncation  error  since  u.  v 
along  with  the  coefficients  e^,  e~ ^  often  vary  sharply  between  discreuzaiion 
nodes. 

2.1  Scherfetter  -  Gummel  discretization 

For  any  2  adjacent  nodes  i,j  of  a  discretization  grid,  the  i,j' th  entry  in  the 
Jacobian  matrix  is  taken  to  be; 


aij]  =  %-d*’B(to  -  to) 

U>3 

(4) 

=  %-d~*’B(to  ~  *j) 

a;. 

(5) 

where  dXJ  is  the  distance  between  the  nodes,  c,y  is  a  cross  -  section  through 
which  the  current  is  assumed  to  flow,  and  B(y )  is  Eernulli’s  function: 

This  scheme  can  be  derived  as  an  exact  solution,  assuming  constant  current 
flowing  parallel  to  the  edge  connecting  the  nodes  i,  i.  and  piecewise  linear 
electrostatic  potential  [5]  [2]. 

2.2  Stream  Functions 

Assume  for  simplicity  R  =  0.  This  assumption  will  be  dropped  later.  If  we 
denote  the  electron  and  hole  currents  by: 


then  (  2)(  3)  can  be  rewritten  as  carrier  conservation  statements: 


<1 

1! 

o 

(9) 

<1 

II 

o 

(10) 

A  solenoidal  vector  field  can  be  expressed  as  a  curl  of  a 
[7],  hence  there  exist  vectors  0n,0p  such  that: 

’stream^  function 

Jn-^xOn 

(11) 

Jp  -  V  x  0P 

(12) 

In  R2  the  vectors  6n ,  9p  point  to  the  third  direction,  and 
be  regarded  as  scalar  quantities  0n,0v. 

consequently  can 

Dividing  (  7).(  8)  by  e^,e  ^  respectively,  and  applying  the  curl  operator 

result  in: 

V  x  e~'JjJn  =  V  x  Vu  =  0 

(13) 

V  X  e'-’Jp  =  V  X  Vt>  =  0 

(14) 

Substituting  (  11),  (  12)  we  obtain: 

V  x  e"*V  X  0n  =  0 

(15) 

V  x  e*V  x  Op  =  0 

(16) 

In  R2  these  reduce  to: 

V  -e^VOn  =  0 

(17) 

V  •  e^Vflp  =  0 

(18) 

where  0n,0p  should  be  regarded  as  the  scalar  magnitude  of  the  corresponding 
vectors. 

Comparison  of  (  2),(  3)  vs.  (  17),  (  18)  reveals  the  following  result: 


Theorem  1  In  R2  and  without  recombination,  the  equation  for  u  is  equiva¬ 
lent  to  the  equation  for  6p,  and  the  equation  for  v  is  equivalent  to  the  equation 
for  0n . 

At  this  point,  if  we  are  determined  to  justify  the  empirical  observation  that  9n 
and  0p  are  smoother  than  u,v  then  the  last  resort  is  the  boundary  condition 
configuration. 

2.3  Boundary  Conditions 

The  nature  of  the  boundary  conditions  for  0n,  0p  is  revealed  via  the  following 
lemma: 

Lemma  1  The  equipotential  curves  of  0n  are  tangent  to  Jn ,  and  normal  to 
the  equipotential  curves  of  u.  The  same  statement  applies  to  0p ,  Jp.  v. 

Proof:  V#n  _L  Vu  since 

V0n  •  Vu  =  V0n  ■  e~'i  Jn  =  V6>n  •  e“*V  x  0n  =  0  (19) 

The  curves  Gn  =  const,  are  normal  to  V0n,  hence  parallel  to  Vu  which  is 
Jn/e'i'  The  curves  on  which  u  =  const,  are  normal  to  Vu  hence  also  to  the 
curves  9n  =  const.  The  proof  for  0p,Jp,v  is  identical.  □ 

Corollary  1  .  6n  is  constant  on  contiguous  segments  ofT jv,  and  d6n/di /  = 
0  on  To-  Similar  statements  apply  to  6p. 

The  resulting  problem  for  0n,0p  is  still  under  determined.  For  each  segment 
of  T.v  the  actual  constant  value  of  0n,0p  has  to  be  found.  The  necessary  ad¬ 
ditional  constants  are  easily  derived  from  the  Dirichlet  boundary  conditions 
for  u,  v.  Let  the  Dirichlet  boundary  be: 

m 

rD  =  U  rb  (20) 

»=i 

Where  each  Tp  is  connected.  Integrating  (  7),  (  8)  we  obtain: 

u(xi)  =  «(x0)  +  Jc  (21) 

u(i1)  =  u(z0)+  [  JP(s)e'l’{3)ds  (22) 

Jc 


$ 
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where  C  is  any  path  leading  from  Zo  to  xj.  Choosing  arbitrary  curves  from 
rj,  to  the  m  -  1  other  segments,  equations  of  the  type  (  21), (  22)  give  m  —  1 
constraints: 


ut-u1=  [  (V  X  ffn)e~*Md3 

JC , 

(23) 

Vi  -vi=  /  (V  x  (i^e^ds 

(24) 

jC, 


There  remains  an  arbitrary  constant  in  the  problem,  which  is  physically 
immaterial.  When  m  =  1  then  8n,8p  are  both  constants,  and  there  is  no 
current  flowing  through  the  device. 

We  can  now  attempt  to  »olve  the  problem  by  the  following,  loosely  defined 
al  gorithm: 

Algorithm  1. 

1.  Discretize  (  17), (  IS). 

2.  Discretize  the  m-1  constraints  (  23), (  24) 

3.  Solve  the  resulting  (under-determined)  linear  system. 

4.  Integrate  lor  the  values  of  u,v  using  (  21),  (  22). 

Mock  [6]  showed  that  the  Scharfetter  -  Gummel  discretization  scheme  is 
equivalent  to  a  finite  element  version  of  this  algorithm,  using  quadralateral 
piecewise  constant  elements  for  8n,8p. 

2.4  Nonzero  Recombination 

If  V  •  J  =  R  ^  0  the  J  is  not  solenoida.l,  and  cannot  be  represented  as  a  curl 
of  a  vector  potential.  However  if  we  definethe  2  current  components  to  be: 

j\{x,y)  =  Mx,y)  +  J  R(t,y)dt  (25) 

h{x,y)  =  Mx^y)  (26) 

then  J  is  a  solenoidal  vector  and  the  preceeding  derivation  can  be  applied 
to  it. 
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3  examples 

з. 1  A  Capacitor 

A  semiconductor  capacitor  is  a  simple  ID  device  insulated  at  one  end.  As¬ 
sume  a  p  -  type  silicon  piece  insulated  on  the  left.  A  contact  attached  to  the 
insulator  is  called  the  ‘gate  contact’,  and  a  contact  at  the  other  end  of  the 
device  is  called  the  ‘source  contact’.  When  a  positive  bias  is  applied  between 
gate  and  source,  holes  are  repelled  from  the  insulator  and  a  depletion  layer 
of  unbalanced  acceptor  ions  is  formed  next  to  the  gate.  At  a  higher  voltage 
there  appears  a  thin  layer  of  electrons  near  the  gate.  This  layer  is  called 
the  ‘inversion  layer’.  Being  extremely  thin  and  steep,  the  inversion  layer  is 
very  hard  to  resolve  by  discretization.  While  being  a  rather  simple  device, 
capacitors  are  often  imbedded  in  more  complex  semiconductor  devices  (par¬ 
ticularly  MOS  devices),  and  the  difficulties  mentioned  above  are  inherited 
by  their  simulations. 

It  should  be  mentioned  here  that  while  n,p  have  a  severe  boundary  layer 
at  the  gate,  u,v  are  actually  constant  for  a  capacitor.  Nevertheless,  the 
conservation  equations  for  either  the  n,p  variables,  or  for  u,u  are  highly 
ill-conditioned  [l]. 

The  solution  of  such  a  capacitor  by  Pisces  [4]  requires  an  extreme  mesh 
refinement  at  the  inversion  layer.  For  gate  to  source  bias  which  exceeds  C.b 
volts,  the  Jacobian  matrices  ocurring  in  the  Newton  process  become  too  ill 
conditioned,  and  Pisces  ceases  to  converge. 

In  order  to  apply  algorithm  1.  to  the  capacitor  we  regard  it  as  a  long  and 
narrow  2-D  device,  with  carrier  Dirichlet  boundary  conditions  at  the  source, 
and  Neumann  boundary  conditions  along  the  remaining  boundary.  While 

и,  v  are  separated  from  their  boundary  data  by  regions  of  extreme  coefficient 
variation,  9n,0p  are  always  near  their  boundary  data.  The  equations  for  the 
stream  -  functions  are  hence  very  well  conditioned.  The  integrals  (  21), 
(  22)  may  produce  large  sums,  but  this  integration  operation  is  perfectly 
well  conditioned,  and  only  requires  an  appropriate  quadrature  rule. 

It  should  be  stressed  that  despite  the  similarity  between  the  Scharfetter  - 
Gummel  scheme  and  the  stream  -  function  formulation,  the  latter  is  very 
well  conditioned  in  the  capacitor  case,  while  the  former  may  be  highly  ill- 
conditioned. 
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3.2  Diodes 


np-diodes  are  characterized  by  a  doping  function  which  is  nearly  piecewise 
constant,  but  varies  steeply  in  a  narrow  region  called  ‘a  junction’,  tx,  u,  or  n,  p 
normally  have  corresponding  steep  gradients  (or  discontinuities)  about  the 
junction.  Dirichlet  boundary  conditioned  for  the  diode  are  specified  at  the 
contacts,  and  Neumann  boundary  conditions  along  the  other  2  edges.  The 
carrier  variables  are  always  separated  from  some  of  their  Dirichlet  data  by 
the  junction.  Using  the  alternative  stream  function  formulation,  we  observe 
that  the  junctions  are  normal  to  the  Dirichlet  boundaries.  6n,6p  are  ‘short- 
circuited'  by  the  boundary  conditions  (  23),  (  24),  and  consequently  vary 
smoothly  through  the  junction.  As  in  the  capacitor  case,  (  21),  (  22)  have  to 
be  integrated  to  obtain  the  carrier  concentrations,  but  using  an  appropriate 
quadrature  this  operation  poses  no  numerical  difficulty. 

3.3  MOS  transistors 

3.4  MOS  devices 

A  typical  2-D  model  for  a  MGS  device  is  a  rectangular  silicon  piece,  part  of  its 
upper  edge  is  covered  by  a  thin  oxide  layer  (an  insulator).  An  npn  type  MOS 
transistor  has  2  small  heavily  doped  n  regions  by  the  top  corners,  while  the 
rest  of  the  device  is  lightly  p  doped.  Source  and  drain  contacts  are  attached 
to  the  n  regions,  and  a  gate  contact  is  attached  to  the  insulator.  We  can 
distinguish  2  types  of  substrate  boundary  conditions  (applied  at  the  bottom 
edge  of  tne  device):  The  substrate  may  be  grounded  (Dirichlet  boundary 
condtions)  or  insulated  (Neumann  conditions).  An  insulated  substrate  leads 
to  a  ‘floating  region’,  namely,  a  region  without  a  contact.  Floating  regions 
often  lead  to  extreme  ill  conditioning  of  the  Jacobian  matrix  [l] 

The  abrupt  interfaces  between  the  p  region  and  the  n  regions  cause  sharp 
variation  of  the  carrier  concentrations,  as  described  for  the  diodes.  When 
a  positive  voltage  is  applied  to  the  gate,  a  depletion  layer  and  an  inversion 
layer  form  next  to  it.  This  thin  inversion  layer  is  called  ‘a  channel"  since  it 
allows  for  a  controlled  amount  of  current  to  flow  between  surce  and  drain. 
In  the  case  of  a  grounded  substrate  Tp  is  not  much  longer  then  Tn,  and  con¬ 
sequently  it  is  not  clear  that  9n,0p  have  more  favorable  boundary  conditions 
than  n,p.  In  the  case  of  a  floating  region  however,  the  stream  functions  have 
Dirichlet  boundary  conditions  almost  around  the  entire  device.  Therefore 
the  application  of  algorithm  1  leads  to  a  well  conditioned  system,  and  is 
numerically  superior  to  the  Scherfetter  -  Gummel  discretization. 
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Notice  that  even  though  9n  has  a  Dirichlet  type  boundary  at  the  gate  in¬ 
terface,  the  fact  that  e^’  is  much  larger  near  the  inversion  layer  than  in  the 
rest  of  the  semiconductor  forms  a  separation  layer,  in  which  the  coefficient 
of  (  l7)  is  very  small.  This  explains  the  boundary  layer  for  9n  which  exists 
in  the  channel. 

4  Conclusion 

We  have  shown  that  stream  potentials  for  the  semiconductor  continuity  equa¬ 
tions  obey  an  adjoint  set  of  boundary  conditions.  This  fact  explains  the  ad¬ 
vantage  of  using  harmonic  averages  for  certain  geometrical  configurations.  El 
conditioning  of  the  Jacobian  matrices  which  occurrs  in  the  discretization  of 
capacitors,  and  floating  region  devices  can  be  removed  by  solving  the  stream 
function  equations  directly. 
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